Three-dimensional Brownian diffusion of rod-like macromolecules in the presence of 
randomly distributed spherical obstacles: Molecular dynamics simulation 
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Brownian diffusion of rod-like polymers in the presence of randomly distributed spherical ob- 
stacles is studied using molecular dynamics (MD) simulations. It is observed that dependence of 
the reduced diffusion coefficient of these macromolecules on the available volume fraction can be 
described reasonably by a power law function. Despite the case of obstructed diffusion of flexible 
polymers in which reduced diffusion coefficient has a weak dependence on the polymer length, this 
dependence is noticeably strong in the case of rod-like polymers. Diffusion of these macromolecules 
in the presence of obstacles is observed that is anomalous at short time scales and normal at long 
times. Duration time of the anomalous diffusion regime is found that increases very rapidly with 
increasing both the polymer length and the obstructed volume fraction. Dynamics of diffusion of 
these polymers is observed that crosses over from Rouse to reptation type with increasing the density 
of obstacles. 



I. INTRODUCTION 

The phenomenon of Brownian diffusion of particles in 
the presence of fixed obstacles - obstructed diffusion - is 
very ubiquitous in biology, chemistry, and physics. This 
phenomenon, as a combination of several concepts such 
as obstructed random walk, hydrodynamic interactions 
and topological hindrances, introduces a chahenging the- 
oretical problem. Two simple models of this problem 
namely three dimensional random walk in a regular lat- 
tice with a fraction of its sites obstructed and diffusion of 
a pointlike tracer in an ordered medium of spherical ob- 
stacles have been studied theoretically and using lattice 
Monte Carlo method, respectively Also, existence 

of anomalous diffusion regime in Brownian diffusion of 
particles due to the presence of fixed and mobile obsta- 
cles has been reported '^--51 . 

In the context of obstructed diffusion, the other sub- 
ject is the diffusion of solute molecules in hydrogels for 
which there are numerous suggestions (analytical, simu- 
lation, and empirical) for diffusion coefficient as a func- 
tion of obstacles density. In a review of studies of this 
phenomenon it has been shown that diffusion of solute 
molecules in hydrogels composed of rigid (flexible) poly- 
mers can be described reasonably by obstruction (hydro- 
dynamic) models [6*]. Also, the effect of the structure of 
obstacles medium on diffusion of fluid molecules has been 
studied recently 0- 

The other issue in this category is the Brownian dif- 
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fusion of polymeric macromolecules in a medium of ob- 
stacles which has different branches depending on the 
polymer chains properties and the strength of induced 
confinement. In diffusion of polymer chains in a ran- 
dom media, such as exclusion chromatography, the poly- 
mers enter in cavities of typical size comparable to or 
larger than their size. In transition between these cavi- 
ties the chains have to overcome entropic barriers which 
are set up due to reduction of their possible configura- 
tions |8l-[ll|. In a semidilute solution of polymer chains 
or when a polymer chain diffuses in a dense medium of 
fixed obstacles, typical size of cavities is smaller than the 
polymer size and both topological hindrances and en- 
tropic effects are of important role p^ - [l^ . Ratio of the 
strength of these two effects depends on the flexibility of 
the chains. The extreme case of a stiff polymer chain is a 
rod-like polymer for which internal entropy is negligible 
and topological interactions are dominant in its diffusion 
in a crowded environment. 

Despite pure translational Brownian motion of small 
isotropic particles, Brownian motion of anisotropic 
macromolecules such as rod-like polymers contains both 
translational and rotational parts. Coupling of trans- 
lational and rotational motions of such macromolecules 
makes understanding and visualization of their Brownian 
motion noticeably difficult [l^ . Anisotropy affects Brow- 
nian motion of a macromolecule more dramatically in 
the presence of other macromolecules or fixed obstacles. 
In this case, both hydrodynamic friction coefficient and 
confinement effects due to the presence of other macro- 
molecules or obstacles are anisotropic. A rod-like macro- 
molecule experiences the least (the most) collisions with 
other macromolecules or obstacles when it moves along 
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(perpendicular to) its axis. 

Numerous interesting phenomena are known to arise 
from special shape of rod-like macromolecules. Nematic 
ordering in dense solution of these macromolecules due 
to only excluded volume interaction (l6l [l7j . a variety 
of their orientational orderings resulted from the com- 
petition of entropic effects with long- and short-range 
interactions 18, 19], and their rich dynamic behavior in 
crowded environments [20l - |23| are examples of these phe- 
nomena. 

In this paper we study Brownian diffusion of a spheri- 
cal particle and rod-like macromolecules in the presence 
of randomly distributed spherical obstacles using molecu- 
lar dynamics simulations. We find that dependence of the 
reduced diffusion coefficient of the spherical particle and 
rod-like macromolecules on the available volume fraction 
can be described by a simple power law function. In the 
case of rod-like macromolecules, the exponent of men- 
tioned power law dependence is an increasing and satu- 
rating function of the macromolecule length. Despite the 
case of flexible polymers in which reduced diffusion coef- 
ficient has a weak dependence on the length of polymer 
[l3| , noticeably strong dependence is observed here which 
originates from stiffness of rod-like macromolecules. In 
obstructed diffusion of these polymers, an anomalous dif- 
fusion regime is observed which its duration time is a 
rapidly increasing function of both the polymer length 
and the obstacles density. With increasing the density of 
obstacles, long time diffusion of these rod-like polymers 
crosses over from Rouse to reptation dynamics. 

The rest of the paper is organized as follows. The 
model and the simulation method are described in Sec. 
im The results are presented in Sec. IIIIl Conclusions 
and a short discussion are presented in Sec. IIVI 



II. THE MODEL AND THE SIMULATION 
METHOD 

In our simulations which are performed with the MD 
simulation package ESPResSo ^24,], a polymer is mod- 
eled as a bead-spring chain of length N {N spherical 
monomers of diameter a). Successive monomers of the 
polymer chain are bonded to each other by a FENE (fi- 
nite extensible nonlinear elastic) potential psj . 
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with bond strength khond ~ SOs/a^ and maximum bond 
length Rq — l.ba. Bending elasticity of the chain is mod- 
eled by a bond angle potential. 



Ubendir) = fc6end(l 



(2) 



in which is the angle between two successive bond vec- 
tors and kbend is the bending energy of the chains. The 
value of the chain persistence length relative to its con- 
tour length, Y', which is a measure of the flexibility of 



the chain depends on the value of fcbend as Ip = -j^a. 
To model a rod-like polymer we consider Ip ~ lOOLc 
in our simulations. We also model the spherical obsta- 
cles by fixed spheres of diameter a which are randomly 
distributed inside the simulation box with their overlap- 
ping permitted. Excluded volume interaction between 
monomers and the obstacles is modeled by a shifted 
Lennard- Jones potential. 
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in which e and a are the usual Lennard- Jones parame- 
ters and the cutoff radius Average dis- 
tance between neighboring monomers which results from 
combination of FENE and Lennard- Jones interactions 
is 6 ~ 0.97ct. Simulation box is cubic with periodic 
boundary conditions and the temperature is kept fixed 
at ksT — l.Oe using a Langevin thermostat. Simulation 
box length is i = 20(7 and L = 25(T in simulations of 
polymers with iV < 17 and the polymer with iV = 21, re- 
spectively. MD time step in our simulations is r = O.OItq 



in which tq = \J time scale and m is the 

mass of the monomers. Dynamics of the monomers is as- 
sumed to be overdamped by using the friction coefficient 
r = ^ for each monomer in the Langevin thermostat. 
In fact, we model the solvent implicitly by the stochas- 
tic and dissipation terms in Langevin equation instead 
of adding solvent molecules explicitly to the system. Ac- 
cordingly, hydrodynamic interactions are not considered 
here and dynamics of polymer chains in the absence of ob- 
stacles is of Rouse type. In each simulation we first ran- 
domly distribute the obstacles inside the simulation box 
in a way that their overlapping is permitted. Then we 
put the rod-like polymer in a random position with a ran- 
dom orientation inside the box and start the simulation. 
After equilibration of the system which is needed because 
of used potentials and fixed temperature, we probe center 
of mass position, orientation, and velocity of the polymer 
for a long time {Nt time steps) to calculate correlation 
functions. From recorded position of a diffusing particle 
(position of the center of mass in the case of rod-like poly- 
mers), mean square displacement (MSD) at time t = rifT 
is calculated as 
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in which f(t) is the position vector of the particle at time 
t. Diffusion coefficient for each realization of obstacles is 
obtained from time dependence of MSD. For fixed values 
of the polymer length and obstacles density we repeat 
our simulations with 10 different realizations of obstacles 
to average quantities and obtain error bars. 
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FIG. 1. (Color online) Reduced diffusion coefficient of a spher- 
ical particle in the presence of randomly distributed spherical 
obstacles as a function of obstructed volume fraction. Some 
suggested functions for dependence of diffusion coefficient on 
obstructed volume fraction are also plotted. Simulation data 
fits well with a simple power law function of available volume 
fraction with exponent b = 0.86. Inset: A sample simulation 
data for particle MSD versus time in logarithmic scales. Bal- 
listic and diffusion regimes can be seen (slopes of dashed and 
solid lines are 2.0 and 1.0, respectively). 
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FIG. 2. (Color online) Reduced diffusion coefficient of the 
center of mass of rod-like polymers of different length versus 
obstructed volume fraction, Data for spherical particle, 
TV = 1, is shown again for comparison. Power low function 
t)cm ~ (1 — (/>)'' can be fitted to simulation data for all values of 
reasonably (dashed lines). Mackie-Meares function is also 
plotted by solid line (see Sec. IIIIBl in the text). Inset: The 
exponent h versus polymer length, TV, that is an increasing 
semi-linear function showing saturation for N > 17. 



III. RESULTS 

To study obstructed diffusion of a particle or a macro- 
molecule, we first obtain its diffusion coefficient in the 
absence of obstacles (-Do) at given temperature, T, and 
friction coefficient, F, which is in consistence with Ein- 
stein relation, Dq = Then we calculate diffusion 
coefficient, D, in the presence of fixed spherical obstacles 
of given density. Obstructed volume fraction, 0, is cal- 
culated using ESPResSo package (a mesh based cluster 
algorithm [2J|). This calculation is checked and justified 
with examples for which obstructed volume fraction can 
be calculated exactly. Reduced diffusion coefficient for 
each value of 6 is defined as D = . 

A. Obstructed diffusion of a spherical particle 

From MD simulations of the Brownian diffusion of a 
spherical particle in the presence of fixed spherical obsta- 
cles, we calculate diffusion coefficient at different values 
of obstacles density. Diffusion coefficient is calculated 
from the plot of the particle MSD versus time when its 
dynamics becomes of diffusion type at long enough times 
(see the inset of Fig. [ij . Reduced diffusion coefficient of 
the spherical particle versus obstructed volume fraction, 
averaged over 10 realizations of the system for each value 
of (p is shown in Fig. [1] Examples of functions suggested 
for dependence of D such as D — {1 — cj))^/'^ form Ref. 

D = l/Xl -f- 0/2) from Ref. [2] and D = {I - (t)f/^ 
form Ref. 12 Tj are also shown in this figure. As it can be 
seen, the best fitted function is a power law function of 



available volume, 1 — 0, with an exponent h = 0.86. The 
next closer function is the prediction of available volume 
theory [27| . Our data points are considerably far from 
that of Ref. at high values of 0. As we checked by 
running the simulation with the same configuration of the 
obstacles as of Ref. [3] , this difference doesn't originate 
from obstacles configuration. This difference may be be- 
cause of the fact that studying the diffusion phenomena 
in the presence of obstacles by lattice and continuous- 
space methods gives different results [7|]. Note that the 
difference between our results and those of Ref. [l^ orig- 
inates from the fact that we have not used explicit fluid 
particles in our simulations. Dq in our work is constant 
and does not depend on fluid density or free volume frac- 
tion. In the case of Ref. ^21] however, Dq depends on 
the fluid density (Eq. 6 of Ref. [13). 

B. Brownian diffusion of rod-like macromolecule in 
the presence of spherical obstacles 

In simulations of the Brownian motion of rod-like poly- 
mers in the presence of spherical obstacles we use Rouse 
dynamics. Accordingly, hydrodynamic interactions are 
not considered and in the absence of obstacles friction 
coefficient is not direction dependent. Anisotropy in dy- 
namics of these macromolecules originate solely from ob- 
struction effects of the obstacles. After equilibration of 
the system for each value of obstructed volume fraction, 
0, we record the position of the polymer center of mass 
and its orientation (a unit vector parallel to the poly- 
mer axis) at all simulation time steps. From trajectory 
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FIG. 3. (Color online) Log-log plot of time-dependent cen- 
ter of mass MSD for a polymer of length A^ = llat0 = O 
and = 0.52. As it can be seen, despite the = case an 
anomalous diffusion regime at intermediate times is obviously 
seen in the case of (/> = 0.52. In time interval between bal- 
listic and diffusion regimes (approximately from i = 0.1 to 
t = V) velocity autocorrelation function is not yet vanished. 
Duration time of anomalous diffusion regime is observed that 
increases with increasing the polymer length and the value of 
Inset: Duration time of anomalous diffusion regime ver- 
sus the polymer length at three different obstructed volume 
fractions. Simulation of the polymer with A*' = 21 monomers 
at (/) = 0.7 takes very long time and we had not a reasonable 
data for these parameters. 



of the center of mass and time dependence of polymer 
orientation we calculate diffusion coefficient and probe 
orientational relaxation of the polymer. Reduced diffu- 
sion coefficient of the center of mass versus (/> for four 
polymers of different length are shown in Fig. [51 These 
diffusion coefficients are obtained from time dependent 
MSD of the polymers center of mass, two samples of 
which are shown in Fig. [31 Considering that functional- 
ity of dependence of the diffusion coefficient on the ob- 
structed volume fraction is an important question in the 
subject of obstructed diffusion, we search this function- 
ality here. Our results show that this dependence in the 
case of rod-like polymers as well as spherical particle can 
be reasonably described by a power law function of avail- 
able volume fraction, (1 — </>)'', in which the value of the 
exponent h depends on the length of the polymer. It has 
akeady been suggested that obstructed diffusion coeffi- 
cient can be written as a function of a single parameter, 
namely the available volume fraction 21\. As it is shown 
in the inset of Fig. [51 h grows linearly with increasing the 
polymer length, iV, and becomes saturated. The function 
for dependence of D suggested by Mackie and Meares 

a, 



1 



is also plotted as solid line in Fig. [2l This function is 
quite far from spherical particle and short polymers data. 
It is close to the data of long polymers of lengths = 17 
and iV = 21. Noticeable dependence of Dcm on polymer 
length which can be seen in this figure is quite different 
from the case of flexible polymers which is reported in ref. 
[l3l |. In a crowded environment that diffusion dynamics 
of individual monomers is slow relative to the case with 
no obstacles, diffusion coefficient of a flexible polymer is 
also diminished by the same factor as monomers diffusion 
coefficient. In this case, the polymer reduced diffusion co- 
efficient becomes independent of its length. In the case of 
rod-like polymers however, the monomers have to obey 
whole polymer dynamics because of its rigidity. In this 
case longer polymers experience stronger obstruction be- 
cause of stronger spatial hindrance and the polymer dif- 
fusion coefficient considerably depends on the polymer 
length. Dependence of flexible polymers diffusion coef- 
ficient in disordered porous material on matrix volume 
fraction when numerous polymers diffuse simultaneously 



has been already reported [14 



(5) 



In the plot of center of mass MSD versus time in log- 
arithmic scales, Fig. [31 at = in addition to ballistic 
and diffusion regimes which correspond to linear parts of 
slopes 2 and 1, respectively, a third semi linear region can 
be seen. In this region velocity autocorrelation function 
is not yet vanished and the diffusion regime is not yet 
started. As it is shown in Fig. [3l another linear regime 
in addition to ballistic and diffusion regimes can be seen 
in which velocity autocorrelation function is vanished but 
the slope of fitted line differs from unity. This region cor- 
responds to anomalous diffusion of the polymer over a 
time interval during its dynamics. Our results show that 
duration of anomalous diffusion of the rod-like polymer 
increases rapidly with increasing its length and the value 
of the obstructed volume fraction, 0. In the inset of Fig. 
[31 time length of anomalous diffusion regim of the poly- 
mer, Tq, versus its length at three different values of 
are shown. These times are roughly obtained from our 
MSD data as the time interval in which data points are 
clearly out of two lines of slopes 2 and 1 and the semi- 
linear part in which velocity autocorrelation function is 
not yet vanished (see main part of Fig. [3]). As it can be 
seen, the value of changes its order of magnitude with 
increasing the strength of obstruction, namely the value 
of and the length of the polymer, IS! . Increasing of the 
duration of anomalous diffusion regime with increasing 
the obstructed volume fraction has also been reported 
already in the context of particle diffusion in quenched 
media (see for example refs. 0, Q)- 

To investigate anisotropy in diffusion dynamics of rod- 
like polymers in strongly obstructed conditions, a sam- 
ple of which is shown in Fig. [31 we calculate diffusion 
coefficients D|| and D_l in directions parallel and per- 
pendicular to the polymer axis, respectively. It is done 
by decomposing displacement of the polymer center of 
mass in each time step into two parts, parallel and per- 
pendicular to its axis. For a rod-like diffusing object it is 
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FIG. 4. (Color online) Sample trajectory of the center of 
mass of a rod- like polymer with A'^ = 11 monomers at two 
obstructed volume fractions (/> = 0.1 and </> = 0.7 during At = 
5 X lO" MD time steps. Trajectory in the case </> — 0.7 is 
more anisotropic and contains linear parts showing that the 
polymer dynamics is of reptation type. 



known that 



Du+2D, 



(6) 



and in the reptation regime where D± ~ 0, this relation 
becomes Z?cm — For a rod-like polymer consisting 

of iV = 17 monomers, and D± are shown as func- 
tions of (j) in Fig. [S] With increasing the value of <j>, 
it can be seen that D± vanishes meaning that the poly- 
mer dynamics becomes of reptation type. It is known in 
polymer physics that despite the Rouse dynamics in 
which Dcm ~ N~^, in the reptation regime, Dcm N^'^. 
In the inset of Fig. [Sj log-log plot of Dcm versus polymer 
length, N, is shown. Crossing over from Rouse to repta- 
tion dynamics with increasing is obviously seen in this 
figure. 

To investigate orientational relaxation of a rod-like 
polymer which is related to duration time of their anoma- 
lous diffusion regime, we calculate autocorrelation func- 
tion of the unit vector corresponding to its direction, u{t), 
defined as 



u{t).u{0) 



1 



Nt-nt 



Nt-nt 



u{{ns + nt)T).u{nsT). (7) 



Relaxation time, r^, is calculated by fitting the expo- 
nential function, Aexp(— t/r,.), to the simulation data as 
is shown in the inset of Fig. [5] for a polymer of length 
Af = llat0 = O.3. Also, following Ref. 15|, we probe re- 
laxation of diffusion coefficients along and perpendicular 
to the polymer initial direction from their initial values to 
the value of Dcm over the time. At small times, diffusion 
coefficient along the polymer is larger than that in per- 
pendicular direction. As time goes on and the polymer 



FIG. 5. (Color online) Diffusion coefficient of a rod-like poly- 
mer of length TV = 17 along (fn) and perpendicular to {D±^ 
its axis versus (j). As it can be seen, D[| and D±_ are respec- 
tively increasing and decreasing functions of 4> showing that 
the polymer dynamics crosses over to the reptation type. In- 
set: Log-log plot of Dcm versus polymer length, A^, which 
shows that with increasing the fraction of obstructed volume, 
the exponent a in Dcm ~ A'^"" changes from 1 (Rouse dy- 
namics) to ~ 2 (reptation dynamics). 



forgets its initial direction, diffusion coefficients parallel 
and perpendicular to the initial direction become of the 
same value (see Fig. [B]). As the figure shows, relaxation 
time obtained from both methods are approximately the 
same. According to our results, for a polymer of given 
length, the value of the orientational relaxation time is 
comparable to the time length of its anomalous diffu- 
sion regime. It shows that the existence of anomalous 
diffusion regime comes from anisotropy in the shape of 
the polymer. A rod-like polymer experiences ordinary 
diffusion regime only at time scales larger than its ori- 
entational relaxation time. For a long polymer at high 
values of obstructed volume fraction this relaxation time 
could be very long and may be longer than simulation or 
experiment time. 



IV. CONCLUSIONS AND DISCUSSION 

In conclusion, MD simulation study of the Brownian 
diffusion of a particle and rod-like polymers in the pres- 
ence of randomly distributed spherical obstacles showed 
that: 1- Obstructed volume fraction dependence of the 
reduced diffusion coefficient of a spherical particle, ob- 
tained from our off lattice MD simulations differs from 
that of lattice based study of obstructed random walk. 

2- Reduced diffusion coefficient of both rod-like polymers 
and spherical particle have a power law dependence on 
the available volume fraction. The exponent of the power 
law function in the case of rod-like polymers is an in- 
creasing and saturating function of the polymer length. 

3- Despite the case of flexible polymers for which a weak 
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FIG. 6. (Color online) Relaxation of the diffusion coefficients 
along and perpendicular to initial direction of the polymer 
to the value of its center of mass diffusion coefficient as time 
elapses. Inset: Autocorrelation function of the unit vector 
along the polymer axis (defined in the text) versus the time 
for the same polymer and fitted function, Aexp{—t/Tr), with 
Tr ~ 225. 



dependence of the reduced diffusion coefficient on the ob- 
structed volume fraction has been reported, this depen- 
dence is noticeably strong in the case of rod-like poly- 
mers. 4- Obstructed diffusion of rod-like polymers is 
observed that is anomalous at short time scales. Time 
length of the anomalous diffusion regime is of the same 
order of the polymer orientational relaxation time and is 
a rapidly increasing function of both the polymer length 
and the obstructed volume fraction. 5- With increasing 
the density of obstacles, diffusion dynamics of rod-like 
polymers is observed that crosses over from Rouse to rep- 
tation type. 

Although MD simulation of rod-like polymers as pre- 



sented here is too time consuming because of using nu- 
merous potentials, such simulations doesn't suffer from 
dependence of results on the selected lattice and can give 
better results. Specially, with overdamped dynamics of 
diffusing particle and working with constant tempera- 
ture, the obstruction effect of the obstacles is considered 
very well. 

The main difference between dynamics of flexible poly- 
mers and that of semiflexible or rod-like polymers comes 
from anisotropy in the shape of stiff polymers. Clearly 
the effects of anisotropic shape of a macromolecule on 
its dynamics becomes more intense in the presence of 
obstacles. Dynamics of a flexible polymer in the pres- 
ence of obstacles is dominantly affected by obstruction of 
monomers dynamics alone. However, in the case of stiff 
polymers in addition to obstruction effect of the obstacles 
on the monomers, topological hindrance originated from 
polymer stiffness also plays an important role. Strong de- 
pendence of diffusion coefficient of rod-like polymers on 
the obstructed volume fraction relative to flexible poly- 
mers and existence of anomalous diffusion regime in dy- 
namics of these polymers are consequences of above men- 
tioned difference. 

Analytical approach to obstructed diffusion of rod-like 
polymers seems a challenging problem. Taking into ac- 
count the size of the obstacles and non zero diameter of 
the polymer in obtaining dependence of diffusion coeffi- 
cient on obstructed volume fraction is not straight for- 
ward and an easy task. 
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